Exact Algorithm for Sampling the 2D Ising Spin Glass 



o 
o 

(N 

-t— > 
o 

O 

o 

m 

"a 

C/3 



I 

o 
o 



> 

in 
in 

o 

ON 

o 



Creighton K. Thomas and A. Alan Middleton 

Department of Physics, Syracuse University, Syracuse, NY 13244; USA 

A sampling algorithm is presented that generates spin glass configurations of the 2D Edwards- 
Anderson Ising spin glass at finite temperature, with probabilities proportional to their Boltzmann 
weights. Such an algorithm overcomes the slow dynamics of direct simulation and can be used to 
study long-range correlation functions and coarse-grained dynamics. The algorithm uses a corre- 
spondence between spin configurations on a regular lattice and dimer (edge) coverings of a related 
graph: Wilson's algorithm [D. B. Wilson, Proc. 8th Symp. Discrete Algorithms 258, (1997)] for sam- 
pling dimer coverings on a planar lattice is adapted to generate samplings for the dimer problem 
corresponding to both planar and toroidal spin glass samples. This algorithm is recursive: it com- 
putes probabilities for spins along a "separator" that divides the sample in half. Given the spins on 
the separator, sample configurations for the two separated halves are generated by further division 
and assignment. The algorithm is simplified by using Pfaffian elimination, rather than Gaussian 
elimination, for sampling dimer configurations. For ?i spins and given floating point precision, the 
algorithm has an asymptotic run-time of 0(71^^^^); it is found that the required precision scales as 
inverse temperature and grows only slowly with system size. Sample applications and benchmarking 
results are presented for samples of size up to n = 128^, with fixed and periodic boundary conditions. 



I. INTRODUCTION 

Materials with quenched disorder, such as spin glasses, 
can have extremely long relaxation times, so that labora- 
tory samples exhibit non-equilibrium behavior over many 
decades in time scale Spin glass materials exhibit 

"aging" , a slow evolution in the magnetic response, for 
example, and non-equilibrium phenomena such as "reju- 
venation", where changes in the temperature can undo 
the effects of aging. As these phenomena take place over 
time scales much longer than the microscopic time scale 
for individual spins, these effects must be due to the col- 
lective behavior of many spins. As analytical work is very 
difficult in disordered materials 0, 01 , numerical simula- 
tions have been important in building a picture of the 
low-temperaturephase of models of disordered spin sys- 
tems (e.g., BHIl). 

Numerical work using direct local Monte Carlo simu- 
lation of the dynamics and equilibration Q indicate that 
models such as the Edwards- Anderson model [l^l possess 
the long relaxation times that are at least necessary to 
start to explain these behaviors. Given the direct corre- 
spondence between simulation time and "experimental" 
time, though, the same long relaxation times that one 
is seeking to understand make such simulations very dif- 
ficult, even though very long simulation times are used 
i- 

Various alternate approaches and approximations have 
been developed to address the difficulties of direct simu- 
lation. These approaches can be used to determine both 
the equilibrium state and how this state is approached. 
When the primary concern is the understanding of the 
equilibrium state, many studies have sought to find the 
ground state of given samples, as many of the prop- 
erties of the low-temperature phase are believed to be 
given by the properties of the ground state (such as the 
sample-to-sample fluctuations in the ground state en- 
ergy or the length-dependent domain wall free energy) 



[111 . Il2l . Il3l . Il4l . llSj l . This direction of research is based 
on developing faster exact methods and accurate heuris- 
tic methods for finding the spin configuration that min- 
imizes a Hamiltonian with fixed random couplings. The 
search for a ground state configuration is closely con- 
nected with combinatorial optimization methods devel- 
oped in computer science, though finite-dimensional spin 
glasses additionally lend themselves to real-space tech- 
niques inspired by the renormalization group [ll|. Equi- 
librium quantities at finite temperature, such as the par- 
tition function and density of states, can be computed for 
the 2D Ising spin glass. The approach to the ground state 
and non-equilibrium properties can then be studied by 
direct simulation or possibly heuristically by real-space 
blocking of the degrees of freedom [l^ . 

We present here an algorithm that extends these ap- 
proaches to allow for exactly sampling the configurations 
of the disordered Ising model on 2D lattices without the 
use of Markov Chain Monte Carlo (MCMC). For n spins, 
this algorithm takes 0{n^^^) steps and in practice has 
a running time that grows only somewhat faster, i.e., 
somewhat more rapidly than L^, at fixed temperature. 
As lower temperatures T require more precise arithmetic, 
the running time grows roughly as . The algorithm is 
based on Wilson's algorithm for sampling planar dimer 
models We use a mapping of the Ising spin glass 

model to the dimer problem for the decoration of the 
graph dual to the spin lattice [H, [l^ . We take advan- 
tage of the regular structure of the square lattice to sim- 
plify the algorithm and also modify the matrix algebra of 
Wilson's algorithm so that the calculation is both simpler 
and more numerically stable. 

This algorithm for sampling provides an opportunity 
to study many outstanding questions for 2D spin glasses 
in much more detail than possible with MCMC computa- 
tions. For example, the dependence of replica overlaps on 
temperature and sample size can be directly computed. 
Correlation functions are easily found: these can be used 
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to study the decay of correlations at finite temperature 
in both Gaussian and ±J models, which differ in some 
aspects at T = 0. The power law decay of spin-spin cor- 
relations are presumed to behave as up to the cor- 
relation length: how rj depends on model and is related 
to thermodynamic quantities such as the heat capacity 
is still not completely understood poj . 

A. Model 

The Edwards- Anderson (EA) spin glass model is a pro- 
totypical model for disordered materials. The EA spin 
glass model has the Hamiltonian 

T^jiS) = ~^J^3S^Sj, (1) 

where the J7 = {Jij} are sample-dependent couplings. 
For example, the Jij can be chosen independently and 
randomly from a Gaussian distribution or from a bimodal 
distribution Jij = ±1 (the ±J model), with mean zero 
and variance 1 in either case. These couplings connect 
two neighboring spins, located at points i and j in the 
sample. The spins Si are Ising spins, i.e., each Si = ±1. 
We will only be able to exactly sample in the 2D case. 
We will study the square lattice of spins in both the case 
of periodic boundary conditions, where the bottom row 
of spins is connected to the top and the left column to 
the right column, and the case of fixed boundary con- 
ditions, where the spins on the boundary of the square 
sample are fixed. A spin configuration {si} = S £ S is 
an assignment of spin values Si to each of n sites i; there 
are 2" possible spin configurations in the state space S. 
A ground state spin configuration Sqs that minimizes 
the Hamiltonian can be found in polynomial time using a 
minimum- weight perfect matching algorithm, if the edges 
(ij) which connect nearest neighbor sites and the sites 
{i} form a planar graph p^ . At positive temperature 
T = the partition function for a given realization of 
disorder J is Zj = J^s' 6xp[— /3'Wj-(S")] and the proba- 
bility of observing a spin state 5* in a sample defined by 
J is Pj{S) = Zj^ cxp[-/3Hj{S)] in equilibrium. 

B. Exact computation of the partition function 

It has long been known that the partition function of 
the 2D ferromagnetic {Jij = 1) Ising model with no ex- 
ternal magnetic field can be found exactly by computing 
the determinant of a matrix derived from the spin lat- 
tice. One type of construction of this determinant uses 
a sum over sets of closed loops on the spin lattice: these 
loops represent the terms in a high-temperature expan- 
sion of the partition function. The first published con- 
struction of these type of loops is that of Kac and Ward 
PH . who directly count the polygonal loops. A technique 
for constructing the relevant matrix for the determinant 



technique is to map the Ising model onto a dimer cover- 
ing problem on a decorated lattice G [lol . [2^ , where the 
spins in the original lattice are replaced by a subgraph, 
a Kasteleyn or Fisher city (a dimer covering is a set of 
edges in the graph such that every node belongs to ex- 
actly one selected edge). The Kasteleyn matrix K of the 
graph G for the dimer problem describes the connections 
between neighboring nodes. This square matrix, which 
is indexed by a numbering of the nodes of G, has non- 
zero entries at locations that are indexed by the two ends 
of a connection between the nodes. Counting the parti- 
tion function for dimer coverings is equivalent to com- 
puting the Pfaffian of the Kasteleyn matrix, where the 
Pfaffian in this case is a square root of the determinant. 
These Pfaffian techniques have been used for the exact 
solution of the pure Ising model in the thermodynamic 
limit 0, Ull, [13] and, e.g., for computing the density of 
states in finite samples. Beale [l^] rewrote the Pfaffian 
in a form that allows for faster direct computation of 
the partition function in a pure ferromagnetic model. As 
the derivation of the correspondence between the parti- 
tion function of the Ising model and the determinant or 
Pfaffian methods for finite samples does not rely on a 
homogeneous coupling constant Jij, these methods can 
also be applied to spin glass samples in two dimensions. 
This correspondence has thus been used to compute di- 
rectly the partition function (and density of states) for 
disordered samples [13, [1^ . Pfaffian techniques can also 
be used to compute degeneracies and correlation func- 
tions in the ± J- model (where couplings are all of the 
same magnitude, but randomly ferromagnetic or antifer- 
romagnetic between neighboring spins) [26| and has been 
used to study the heat capacity of this same model at low 
temperatures (e.g., see [27|). 



C. Review of configuration sampling 

Being able to compute the partition function (and of- 
ten the density of states as a by-product) is useful in 
computing such quantities as domain wall free energies, 
sample-to-sample fluctuations in the free energy, specific 
heat, and other global quantities. By computing the 
partition function for fixed relative spin co nfig urations, 
one can also calculate correlation functions [2y] . But for 
many purposes, such as faster computation of correla- 
tion functions, the organization of states in a spin glass, 
or for use in a heuristic for studying the dynamics of dis- 
ordered materials [IB] , it is useful to be able to generate 
sample configurations, given a realization of the disorder. 
For sampling the equilibrium behavior of the system, it 
is sufficient to generate such samples with their proper 
Boltzmann probability Pj{S). For nonequilibrium dy- 
namics, such sampling can be used in patchwork dynam- 
ics, which is closely related to the renormalization ap- 
proaches to nonlocal dynamics used in multigrid Monte 
Carlo methods and hierarchical genetic methods [Tl[.[28t. 

Heuristic sampling, where there is no proof of exact- 
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ness, is typically done using the Markov chain Monte 
Carlo (MCMC) method. In MCMC methods, local prob- 
abilistic dynamics that obey detailed balance are used to 
update the spins. At long times, the probability of ob- 
serving a configuration should be the equilibrium proba- 
bility. The equilibration times using this method can be 
prohibitively long, though, especially in glassy systems 
such as the 2D spin glass [2^. Some faster Monte Carlo 
methods have been developed for the 2D spin glass at low 
temperature (30l | , but with any such method there is also 
a question of how to test whether equilibrium is achieved 
with sufficient accuracy. It is of use to have criteria to 
confirm converges of the Markov chain to the equilib- 
rium distribution. Propp and Wilson [3l| proposed a 
technique for generating exact samples with MCMC by 
"coupling from the past" (CFTP). In this framework, it 
is possible to verify that the system has converged from 
all possible initial conditions to a single state, at which 
point it is exactly in equilibrium. This approach often 
makes use of a natural partial ordering of configurations 
that is used to guarantee convergence. For disordered 
models, there is often no such obvious partial ordering 
of the states that ensures convergence of CFTP. Chanal 
and Krauth |32| have nevertheless succeeded in applying 
CFTP to the Ising spin glass using a coarse-grained orga- 
nization of the states: at first, all states are possible; as 
the Markov chain is developed and the number of states 
is reduced by coupling, the constraint on allowed states is 
further coarse-grained, until a single whole sample state 
is left. But the coupling time (time for convergence to 
a single sample) is still of the order of the equilibration 
time, which of course can be very long at low tempera- 
tures. 

Sampling with the exact Boltzmann weights has been 
implemented and applied to the Migdal-Kadanoff (MK) 
lattice, which is not a finite-dimensional lattice, but is 
used to approximately represent finite-dimensional lat- 
tices. As the MK lattice has a hierarchical structure, the 
spin configurations can be summed over successive scales, 
starting from the smallest, to compute the partition func- 
tion and the relative partition functions can be used to 
sample the spins. This was done in Refs. [H and [s^ to 
study chaos and spin overlap on hierarchical lattices. 

Exact sampling of configurations can always be carried 
out in time polynomial in the size of the sample, if the 
partition function may be calculated efficiently. One di- 
rect, but somewhat slow method, is to assign a single spin 
at random and then compute the partition function con- 
ditioned on assignment of individual neighboring spins; 
this requires n = L'^ computations of the partition func- 
tion for 0{n) spins. Such a technique is mentioned as a 
possibility, for example, in Ref . [ssl. As the partition func- 
tion can be computed in 0(n^/^) steps, this would require 
0(n^/^) arithmetical steps. There are other methods for 
carrying out exact sampling, however. 

Exact sampling of ferromagnetic Ising systems (in any 
dimension) may be performed in polynomial time [36|. 
This technique works in the Fortuin-Kastcleyn cluster 



representation and successively removes bonds and spins 
through a reduction technique. A related problem, sam- 
pling configurations of dimer coverings on a planar bipar- 
tite lattice, has an elegant sampling technique [stI. |38| . 
which exactly maps the statistical mechanics on an L x L 
lattice to an (L— 1) x (L— 1) lattice with modified weights 
on the edges. Other techniques for calculating the exact 
partition function of the 2D Ising Spin Glass, such as 
the Y-A technique of Loh and Carlson [3^, are quite 
similar in spirit to the dimer covering algorithm. This 
technique also involves an efficient recursive reduction of 
any planar graph to a smaller graph, but when frustra- 
tion is present the intermediate reduced bond strengths 
can become complex, which complicates possible sam- 
pling techniques. 



D. Overview of algorithm 

We now outline the crucial points for our application. 
In two dimensions, there is a one-to-one correspondence 
between spin configurations of the Ising model with ar- 
bitrary couplings and dimer configurations on a deco- 
rated version of the dual lattice. The individual spin 
and dimer configurations have the same energy, so the 
corresponding configurations have the same Boltzmann 
weights ex-p{—f3E), where Z and E are the partition 
function and configuration energy for either the dimer 
or spin problem. We can therefore generate sample spin 
configurations by sampling among dimer configurations 
and mapping them to the spin representation. Note that 
the traditional method for calculating the partition func- 
tion is a mapping between the primal lattice and a dimer 
model: a dimer configuration, which defines loops in a 
high temperature expansion of the partition function, 
does not directly map onto a unique spin configuration. 
Using the dual lattice, however, allows for such a map. 

Wilson's algorithm may be used to sample dimer con- 
figurations efficiently for any planar lattice, so efficient 
sampling of the Ising model can be carried out on gen- 
eral planar samples. One requirement for Wilson's al- 
gorithm is an efficient method to recursively subdivide 
the lattice; this task is straightforward on a regular lat- 
tice: we subdivide or separate the sample by choosing 
two adjacent rows or columns of spins. The spins on 
these two lines are the separator sites for the spin lattice. 
These separator spins are then assigned by a sequence 
of weighted choices. The weights for the choice of these 
spins are found, in essence, by computing the needed cor- 
relators between each pair of spins situated on these two 
lines. Once the spins on the separator have been chosen 
and fixed, this division and sampling is repeated on finer 
and finer spatial scales, using the solved spins as fixed 
boundary conditions for the subsamplcs. Besides allow- 
ing for recursive assignment of spins on the separators, 
this nested dissection is used to efficiently organize the 
needed sparse matrix computations. 

We have also simplified the algorithm significantly by 
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using Pfaffian elimination, rather than Gaussian eHmina- 
tion. Pfaffian elimination was used by Galluccio, Loebl, 
and Vondrak (25| in computing the partition function, 
but it can also be used to advantage in sampling. We 
use a sparse matrix representation that greatly reduces 
the amount of space and time needed: due to the regu- 
lar nature of the lattice, all of the primitive operations 
can be explicitly precomputcd and then applied to many 
distinct samples of the same size. We find that the num- 
ber of relevant matrix elements (out of the full O(n^) 
potential elements) that are "visited" during the compu- 
tation scales approximately ^ n and that the number of 
operations obeys the expected growth ~ v?/"^ . 

Though the form of the algorithm that we use is based 
upon and parallels Wilson's algorithm, we present the 
method in detail here. We do this in order to review the 
method itself, emphasize the relationship between match- 
ings and the Ising model, present our form of the matrix 
algebra that we use for sampling dimer matchings, and 
describe sampling for non-planar graphs, such as used for 
periodic boundary conditions. 



E. Implementation results 

As one of the primary motivations for the development 
of our algorithm is its potential use in patchwork dynam- 
ics [T6j . we test our algorithm by timing it in this context, 
random patches of a sample with Gaussian bonds, where 
the variance of the couplings Jij is unity and the mean 
coupling Jij = 0. Our code was developed with the possi- 
bility of using different data types as the matrix elements 
in the calculation. Specifically, we test the algorithm us- 
ing double precision numbers, floating point numbers of 
arbitrary precision, and with exact rational Boltzmann 
weights. As the weights in the computation can vary 
over a large range and a Pfaffian elimination technique is 
used to cancel out matrix elements, similar to Gaussian 
elimination, the algorithm can produce unstable results 
using the floating point types, if proper care is not taken. 
The likelihood of an instability increases with increasing 
system size and with lower temperature. In trying to bal- 
ance the stability and accuracy of the sampling against 
the running time, we determine the arithmetical preci- 
sion needed to reliably sample a configuration. Sample 
results for configurations are displayed in Fig. [TJ Details 
of the precision requirements and example running times 
are given in Sec. IIII El 



II. MAPPING THE ISING MODEL TO A 
DIMER MODEL 

In order to sample Ising spin configurations via the 
sampling of dimer configurations, one requires a one-to- 
one correspondence between the Ising spin configurations 
5 on a given lattice and the dimer covering configura- 
tions Al on a related graph G. Such mappings have been 



constructed for application to the more straightforward 
problem of computing the partition function. These map- 
pings link the problem of computing Zj to & weighted 
enumeration of all perfect matchings A4 on G. A single 
perfect matching on a graph G = (V, E) , where V are 
vertices (nodes) and E are edges connecting pairs of ver- 
tices, is a choice of a subset of edges M C E, the match- 
ing or dimer covering, such that every vertex belongs to 
exactly one edge in M (see Fig. [2|). The generally es- 
tablished procedure for constructing a mapping between 
spin configurations and perfect matchings is to identify 
closed loops on some relevant graph, Go, where Go is ei- 
ther the primary grid (the spin lattice) or the dual lattice 
(the lattice of plaquettes). The partition function, origi- 
nally a sum over spin configurations, can be represented 
as a weighted sum over choices of loops in Go- This 
summation over loops can be carried out by summing 
over matchings on a graph G, constructed by replacing 
the nodes of Go with either Kasteleyn or Fisher "cities" 
subgraphs constructed of a few nodes and edges. 
Perfect matchings on this decorated lattice G then have 
the property that an even number of the covered edges 
are incident upon any given city. The edges of a match- 
ing M that connect cities are therefore even at each city; 
contracting the cities back to single points then gives the 
city-connecting dimers that compose the loops in Gq (see 
Fig.©. 

One mapping between spin configurations and sets of 
loops is based on a high temperature expansion of the 
partition function of the Ising model, where Go is the 
spin lattice and the loops, composed of bonds connecting 
nearest-neighbor spins, represent individual terms in the 
expansion of Zj in powers of exp(— /3Jij). The direct 
replacement of each Ising spin with a "city" gives repre- 
sentation of loops by a dimer matching [l^, [22, . The 
weight of dimer co nfig urations can then be summed using 
Pfaffian methods giving, for example, the Kasteleyn 
solution of the Ising model. However, there is no direct 
correspondence between individual sets of loops and spin 
configurations. 

Alternately, a mapping to G can be defined by taking 
Go to be the dual lattice [H, H^l . This mapping, in con- 
trast with the approach of decorating the original lattice, 
allows for direct sampling of Ising spin configurations. 
The loops on the dual graph represent a loop expansion 
in terms of domain walls. The expansion in domain walls, 
if expressed relative to the ground state, would be a low- 
temperature expansion. More generally, the summation 
is over relative domain walls between a reference config- 
uration and any other configuration. A direct correspon- 
dence between spin configurations and dimer configura- 
tions therefore exists as domain walls uniquely define a 
spin configuration, given a reference configuration, up to 
the possibility of a global spin-flip symmetry. 

Let R = {ri} be a reference configuration of Ising spins 
ri = ±1. We emphasize that this choice is completely 
arbitrary: it need not be a ground state. For convenience 
R can be a configuration with all spins up or a previously 



5 




r=o.5 T=o.2 T=om r=o.02 



Figure 1: Results of applying the sampling algorithm to an individual 2D Ising spin glass sample, for temperatures T = 
0.5, 0.2, 0.08, 0.02, for a single Gaussian spin glass sample with fixed boundaries. The images show the variability of the spin 
assignments (top) and of the domain walls (bottom) over a range of temperatures, in a sample with n — 126^ variable spins 
surrounded by a layer of fixed spins. At least 240 samples were generated at each temperature. The gray scale values indicate 
the probability of a given spin being fixed (upper row) or of neighboring spins being fixed relative to each other (lower row). 
For spin assignments, the darkest colors indicate that the spin is equally likely to be up or down, while light colors indicate 
that the spin occurs with a single alignment in nearly all sampled configurations. These alignments result from correlations 
with the fixed boundary spins. For the domain walls displayed in the lower row of images, the lines indicate the probability 
of relative domain walls between two configurations: the darkest lines indicate the bond dual to that domain wall has a 50% 
chance of opposite or equal relative orientations; where there is no line separating two spins (or only a very light one), the two 
spins have a very high probability of a single relative orientation, either aligned or opposite. Specifically, the bond satisfaction 
variance ^J.i,j{l — is plotted along each dual edge, where ^i^j is the frequency of the JijSiSj being positive. Note that as 
T decreases, the frequency of specific droplet excitations, outlined by domain walls, can either increase or decrease, reflecting 
the sensitivity of the configurations to temperature. This can be seen, for example, in two of the regions that are active at 
T = 0.02, the approximately 20 x 20 region in the far upper left and the approximately 30 x 60 region at the center right: 
the spins in the former become more fixed as temperature decreases while the spins in the latter region become more variable 
(darker) when the temperature is decreased from T = 0.08 to T = 0.02. 



sampled configuration. For a given sampling S of the 
spin configuration, S = {s^}, the loops of dual edges that 
separate spins i and j vifith ViSi ^ rjSj define the relative 
domain walls between R and S. (For the ferromagnetic 
Ising model, one usually takes = 1, so that the domain 
walls separate regions where Si = 1 from regions where 
= -L) 

In this reference configuration, for each pair i, j, define 
Rij = I'ifj as the reference satisfaction of bond i, j . Then, 
for this fixed Rij , we can simply rewrite the Hamiltonian 

as 

T^jiS) ~ — ^ Jij {siSj — Rij + Rij) 

(v) 

= Ur+Hg, (2) 

with I-Lr = — JijRij^ the energy of the reference 

configuration and T-Lq = —J2{ij) Jij i^i^j ~ which 



will be rewritten as the Hamiltonian of the corresponding 
dimer model is the energy of the domain walls between 
the configurations R and S. Note that Ti, r is the same for 
all spin configurations, but must be tracked if comparing 
the effects of changing boundary conditions or comparing 
with ground state energies, for example. 

Let the decorated graph G = (V, E) have the vertex 
set V , which has size \V\ = 2N, N being the number 
of dimers in a perfect matching of the vertices, and the 
edge set E = {cgr} where each edge connects two nodes, 
Ggr = (97'') I for some q,r S V. Then, given a set of 
relative domain wall loops, the dimer configuration is 
uniquely defined by selecting dimers that connect cities 
and cross bonds Jij where SiSj ^ Rij, i.e., that overlie 
the domain walls in Go, and the subsequent unique choice 
of matching for dimers internal to the cities. Choosing 
an energy function w{e) for edges in E with w{eqr) = 
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Figure 2: [color online] A depiction of the correspondence be- 
tween domain wall loops for an Ising spin system and dimer 
matchings on the decorated dual lattice G. (a) A spin sys- 
tem with fixed boundary conditions; an up arrow at location 
i indicates Si = +1 and a down arrow indicates Si = —1. The 
dual lattice Go is indicated by the lines connecting the dual 
nodes, (b) A Fisher city replacement. Each dual lattice node 
is expanded to a Fisher city, a set of six nodes composed of 
two linked triangles, to generate the decorated lattice. For 
work on the square lattice, the bond strengths are set to be 
w{eij) — inside the city, and the bond strengths between 
the cities, indicated here by the notation Wd, d = 0,1,2,3, 
are set according to Eq. ((3|. (c) An example dimer covering 
(i.e., perfect matching) M on the decorated graph G. The 
thicker (also red) bonds with circular ends indicate edges in 
M. The domain walls, composed of dimers that connect dis- 
tinct cities, are indicated by dashed lines, (d) When the cities 
are contracted out from G, the loops on Go remain. Given 
this choice of dimer covering M, the spins that are inside the 
domain walls are flipped to create the new sampled configu- 
ration. 
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Figure 3: The node indexing and edge orientations within a 
Fisher city (left) and the corresponding elements of the 6x6 
submatrix of the Kasteleyn matrix K (right). The numbering 
of nodes is shown for the first city listed in the dual lattice; 
subsequent cities have multiples of 6 added to their indices. 
In the case of the square spin lattice (indicated by the outer 
square bonds on the original spin lattice), all non-zero K ele- 
ments, including d, are set to unit magnitude. The labeling of 
the 0^2 and 1^2 edges indicate how the strengths can be 
modified in the case of the triangular lattice: in this case, one 
can set d — exp[— /3w(ed)] to account for the diagonal bond 
Ed perpendicular to the 2 — > 3 edge. The Kasteleyn matrix 
has row a and column b indices a,6 = 0,...,5. 



for bonds in the cities and 



w{e) 



2 J ij Rij 



(3) 



for dual edges e crossing bonds between spins i and j 
gives 



(4) 



as a consistent energy function for matching configura- 
tions in M. The Ising model and matching model can 
therefore be made equivalent, up to a global energy shift 

Because each dimer configuration corresponds to a spin 
configuration with the same energy, picking a sample 
from the dimer model with the correct probability di- 
rectly produces a corresponding spin configuration that 
has the same probability of occurring. We chose to use 
Fisher cities for this work, instead of Kasteleyn cities p^ . 
as they are simpler to sample using Wilson's algorithm 
on a square lattice. Also, by modifying the weights of the 
Fisher cities, we can also very easily change the weights 
to simulate triangular lattices (see Fig. [3]). 



A. Matchings and the Kasteleyn Matrix 

Given the mapping between matchings using dual lat- 
tice cities and spin configurations, we now briefly review 
the correspondence between dimer matchings and Pfaffi- 
ans. Extensive discussion and examples can be found 
in, for example, Refs. [l^, [s^, [4l[. As a mathemati- 
cal object, the Pfaffian Pf(j4) can be defined for general 
2N X 2N antisymmetric square matrices A = {aqr\q, r = 
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by a restricted sum over per- 
mutations P = p{t) of the indices t = 0, 1, . . . , 2N — 1, 

Pf(A)= (-ir*''^a9ir,ag,r,...a,„r„ , (5) 

P ordered 

where cr(P) is the sign of the permutation from the se- 
quence 0, . . . , 2N — 1 to the sequence qi,ri, . . . , gjv, tn 
and the restriction to ordered P is to rearrangements 
where qk < r^, for all 1 < A; < iV, and qi < q2 < ■ ■ ■ < 
qN- We also have that [Pf(A)]2 = det(^). 

It turns out that summing over permutations with 
these two restrictions is exactly the way to sum over 
dimer coverings for a planar graph G, if the matrix 
elements of A are chosen properly. A matrix whose 
Pfaffian is Zm = <^^p[~/3^g(-^)] is the Kaste- 

leyn matrix K. This matrix has entries K(q,r), with 
q,r = . . . ,2N — 1, satisfying \K{q,r)\ = Xq^r, where 
Xq^r = exp{— /3?i;[e((j', r)]}, and w[e{q,r)] is the bond 
strength associated with edge e{q,r). Directions for the 
edges are then chosen so that all loops in G which en- 
close an even number of nodes include an odd number of 
counterclockwise edges [l9| . The matrix entry K{q,r) is 
set to be Xq^r if an edge is oriented from q to r, otherwise 
it is set to be — Xg,,-- This convention ensures that each 
valid dimer configuration has positive net weight. The 
Kasteleyn matrix is thus a weighted version of a directed 
adjacency matrix. Using these conventions and weight 
assignments gives [l^ 

Pf(A-)= X{xe=ZM=Z^^Zj. (6) 

When decorating Gq with cities to create G, the edges 
internal to the cities must be assigned orientations. An 
example of a Fisher city with the correct directionality 
and the corresponding submatrix is shown in Fig. [3] The 
orientation of the connections between the cities are from 
the 4-node in one Fisher city to the 0-node in the city to 
the right and from the 5-node a city to the 1-node in the 
city in the row above. To simplify notation for the rest 
of the paper, we will use Z to indicate Zq. 

Established analytical and numerical techniques can 
be used to compute Pf(_ft') = Z. As these numerical 
techniques require a number of mathematical operations 
polynomial in the size of the lattice, specifically grow- 
ing as ~ n^/^, the thermodynamic properties can be ef- 
ficiently computed. The number of bits needed for ex- 
act computations grows with n, so that computing, for 
example, the exact partition function, written out as a 
polynomial in exp(— /3) of a spin glass sample for the ± J 
model, where Jy = ±1, requires 0{n7^^) primitive fixed- 
word- length operations [25|. 

We extend this correspondence to carry out sampling 
of spin configurations by applying Wilson's algorithm. 
Partial diagonalization of the Kasteleyn matrix gener- 
ates correlation functions for the choice of the dimers in 
the matching representation. These correlations are be- 
tween dimers on a separator of the sample, which divides 



the sample into two nearly equal places. These corre- 
lations functions include the probability of choosing any 
dimer in a matching, so it is straightforward to determine 
whether a single dimer is selected in a random matching. 
The insight developed by Wilson was to update these 
correlations as dimers are chosen: the effects of partial 
assignment are propagated inductively to correlations be- 
tween other dimers, allowing many dimers to be assigned 
without another factorization of the full Kasteleyn ma- 
trix. Once the dimers have been selected on a separator, 
the two pieces are then solved recursively, using their own 
separators. 



III. WILSON'S ALGORITHM 

In this section, we describe our implementation of Wil- 
son's algorithm, as applied and adapted to sampling con- 
figurations of the Ising spin glass. Wilson's algorithm 
samples dimer coverings: we map the Ising problem to 
the dimer sampling problem using the mapping described 
for the dual lattice in Sec. |TT1 Wilson's algorithm uses a 
"nested dissection" [42|, i.e., a recursive subdivision of 
the sample, where each subdivision of n spins is into two 
pieces of similar size separated by a line of vertices of size 
0{^/n), for efficiency. Such a nested dissection was used 
by Galluccio, Loebl, and Vondrak [2^, to compute the 
full expansion of the partition function of the ±J spin 
glass as a polynomial in exp(— /3), using the high tem- 
perature expansion formulation of the partition function. 
This dissection can be phrased using either a dimer de- 
scription, based on a matching of the decorated graph on 
the dual lattice, or using spins. The algorithm is neces- 
sarily implemented in terms of the former language, but 
for clarity, it is also convenient to describe it using the 
latter language, i.e., based on the spins on the original 
lattice. 

Consider a subsample U of Ising spins {si\i e U}, pos- 
sibly with external fields at the boundary (corresponding 
to fixed spins bordering U ; this graph is still planar). To 
divide this sample into two independent samples, U' and 
U" , a set D of spins is chosen as a spin separator, so that 

U = U'UDUU" (7) 

and no bonds connect spins in U' to spins in U" . We 
choose this spin separator to be composed of two parallel 
lines of spins, so that a line of nodes in the dual lattice 
is contained between the two lines of spins. 

It turns out that Wilson's approach provides an ef- 
ficient way to assign spin values along this separator, 
such that the spins are selected with the correct prob- 
abilities. That is, let such a spin assignment on D be 
Sd = {sk = ±l|fc G D}. The spin at site i for a choice 
5^1 is also written as Soii)- One requires that the proba- 
bility that the algorithm will generate a particular choice 
Sd is just equal to the probability Pj{S\Sd) that the 
properly weighted choice of all spins will yield that par- 
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ticular assignment of spins on the separator _D, i.e., that 



p{Sd) = [zm- 



E 



exp[-/37^(5[/)] 



Su\Su{i)=So{i)yieD 



(8) 

with Sjj being a particular configuration of the spins 
in U, the sum indexing ah possible spin assignments 
consistent with the choice Sd, and with Z{U) = 
exp[— f3H{Su)] the partition function for U. The 
remarkable property of the algorithm to make such a se- 
lection implies that this procedure may then be repeated 
on the remaining unassigned subsystems U' and U" in- 
dependently of one another. 

We can select the assignment for the spins in D by 
sampling from the dimer assignments for all the nodes in 
A, where A is the set of all nodes in G that lie inside of D 
and the connecting edges contained within D. This set of 
nodes A is what is referred to as the separator in Wilson's 
work on an algorithm for random dimcr assignments. 

In order to outline of our version of the algorithm for 
assigning matchings in A, one needs the notion of Pfaffian 
elimination [2^. Let if be a 2N x 2N skew-symmetric 
matrix, i.e., K{q,r) = —K(r,q) for < q,r < 2n ~ 
1. A cross operation between q and j is the addition of 
a multiple of row q to row r and the same multiple of 
column q to column r. If this multiple is given by the 
factor a, the cross operation on K can be written as 



K L{a, q, r)KL'^{a, q, r) . 



(9) 



where L{a, q, r) is the lower triangular matrix / + aSq^r- 
The matrix 6q^r bas all entries zero except for a unit en- 
try in row q and column r. It turns out that the value 
of Pi{K) is unchanged by cross operations, as L has unit 
determinant and Vi{BKB'^) = det(B) Pf (if ) for general 
B [4^ . Pfaffian elimination is the application of multiple 
cross operations to simplify the matrix. This factoriza- 
tion via Pfaffian elimination has the goal of making the 
Pfaffian trivial to compute; the simplest form of a skew- 
symmetric matrix has non-zero values only in the even 
row superdiagonal elements. 



(10) 



1=0 



where is just the matrix that is non-zero except for 
the (2f , 2£ -I- l)'st entry, which is set to 1, and the (2£ + 
1, 2£)'st entry, which is set to —1. In Pfaffian elimination, 
then, the v factors am and the cross operation locations 



are all chosen sequentially so that 
Y = LKL^ , 



(11) 



with L = Y\^m,=i^i'^m,q_rrnTm) and Y is of the form in 
Eq. pU]) . The needed choices of a,„, g„i, and r,„ are 
discussed in more detail in Sec. IIIIBI 

As the factorization of K given by Pfaffian elimination 
leaves the Pfaffian invariant 

Pf (y) = Vi{LKL^) = det(i) Pf (A') = Pf (A') , (12) 



the Pfaffian of the Kasteleyn matrix, and hence the par- 
tition function, can be directly found by multiplying the 
even superdiagonal entries of Y . 

This elimination procedure resembles the application 
of Gaussian elimination to compute the LU factorization 
of a matrix A, with A ~ LU where L is lower triangu- 
lar with unit elements on the diagonal and U is upper 
triangular. The product of the diagonal elements of U 
gives det(A); here Pf(Ar) is the product of the even row 
superdiagonal elements of Y. Factorization via Pfaffian 
elimination maintains the skew symmetry of the par- 
tially factorized l\^L{am,qm,r„i) Kfl^^L'^ {an, qn,rn) 
at each stage. Wilson presented his sampling algorithm 
using Gaussian elimination; we find that Pfaffian elim- 
ination both clarifies the algorithm and makes the pro- 
gramming of the algorithm more direct. A version of 
the algorithm that we implemented using Gaussian elim- 
ination was much less stable numerically than the one 
implemented using Pfaffian elimination. 

The factorization of K given by Pfaffian elimination 
allows the inverse of AT to be quickly computed. It is 
clear from Eq. ([T^ that 



(13) 



where, given the simple form of Y , the inverse of Y is 
easily found: 



Y- 



1=0 ^ 



"2 



(14) 



When the matrix K is created, the indexing of the 
nodes in G is chosen according to a nested dissection of 
the graph G that maintains the grouping of the Kasteleyn 
cities. This ordering reduces the amount of work needed 
to carry out the Pfaffian elimination and is chosen so 
that the elements of the separator at each level of the 
dissection are in a block at the lower right part of the 
submatrix organized by that separator. An example of 
this ordering, given by the nested dissection, is shown in 
Fig.Hl 

The core of the dimer assignment procedure is based 
on the relationship between restricted partition functions 
and the Pfaffian of submatrices of the Kasteleyn ma- 
trix. Consider two partition functions, the entire par- 
tition function Z = Pf(A') and the restricted parti- 
tion function Zp, which is sum of weights OeeGVp ^(^) 
restricted to matchings that include the fixed partial 
matching p ~ {gi, ri, . . . , 5^, r^}, with matched edges 
(gi, ri ),..., (gfc, Tfc). A listing of the terms that con- 
tribute to Zp can be found by removing all nodes in p 
from the graph G and computing the Pfaffian of ifp, the 
Kasteleyn matrix for G\p. To find Zp, the weights x of 
the removed edges must then be included, giving 



Zp = Pf (Ap) Y[ x{e) 



(15) 



The weights x{e) are uniform in Wilson's description, 
though he noted the possibility of variable weights. The 
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(a) (b) (c) 

Figure 4: [color online] An example of the nested dissection and the Kasteleyn matrix K for a 6 x 6 spin lattice sample 
surrounded by an outer layer of fixed spins, (a) The set of 8 x 8 Ising spins sit on the sites of the light gray lattice of bonds 
of strength Jij , where the diagonal bonds are indicated for the case of a triangular lattice. The graph G on which the dimer 
sampling is computed is shown by the darker lines and circular nodes. The gray bands indicate the nested dissection used for 
these nodes: the lighter gray region contains the dimer separator A C G and is bordered by the middle two rows of spins, the 
spin separator D. The darker and medium bands indicate, in order, the subsequent subdivisions of the sample, (b) A display 
of the non-zero elements of the Kasteleyn matrix K, for a left-to-right and top-down ordering of the Kasteleyn cities. The 
nonzero elements of the 294 x 294 Kasteleyn matrix for this are shown as black dots. The edges internal to the Kasteleyn cities 
are closest to the diagonal: further non-zero elements represent connections between the cities, (c) The permuted matrix K, 
where the cities are indexed according to a nested dissection. The gray regions in K include connections contained within each 
of the separators of the nested dissection, with the same shades as in (a) , and between the separators and other nodes. The 
procedure of Pfaffian elimination can at most affect elements within the gray regions and also the values near the diagonal, for 
the nodes not contained in the gray regions in (a). Spin values are assigned to D by examining the part of indexed by 
the nodes of A, i.e., the lower right square submatrix contained within the light gray region. 



probability P{p) of choosing the edge set p is therefore 



(16) 



Given that one has already chosen an edge set p that 
partially covers a graph, the conditional probability 
P{p, u\p) of edge u being in a complete matching that 
includes p is 



Pip, u I P) 



(17) 

Fundamental relations between determinants and in- 
verse matrices are used in Wilson's algorithm to speed 
up the computation of Kp: we directly adapt these re- 
lations for Pfaffian factorization. Let A be a 2m x 2m 
skew- symmetric matrix, and < £ < 2m be an even in- 
teger, and p = {<i,t2, ... be a subset of indices for 
the rows (columns) of A. We will use the notation that 
Ap = At^^t2,...,ti denotes the (2m — I) y. (2m — I) skew- 
symmetric matrix given by removing from A all rows and 
columns with indices in the set (ii, . . . , i^). The notation 
[^]ti,...,t« will denote the £ x £ matrix resulting instead 
from keeping just those rows and columns and eliminat- 
ing the rest of the matrix. Using this notation, and the 
result that det(A) = [Pf(A)]^, Jacobi's theorem (or di- 
rectly using the definition of the Pfaffian to show that 



element i,j of A"! is (-1)'+^ Pf(Ajj)/ Pf(A)) implies 
that H 



Pf(Aii,....if ) 



= ±Pf([A-l], 



(18) 



where the sign depends only on the choice of the indices 
ii, . . . , if. 

Eqs. P?|) and (fT5|) thus allow one to compute the 
probability of matching (gi, ri ),..., (g^, rfc). using the 
Pfaffian of the inverse of the Kasteleyn matrix where 
the same rows and columns kept. The Pfaffian factor- 
ization of this latter matrix, [X^^l , can be 

L J i3i,ri,...,(}fc,rfc ' 

updated incrementally as successive choices of matched 
edges are made. This update allows for the progressive 
computation of the probabilities P{p,u \p) = x{u)zf^[u)^ 
where the updated factorization directly gives the value 

z,{u) = Pi([K-X;) /Pi([K-X). 

Our adaptation of Wilson's algorithm can now be sum- 
marized in outline form: 

1. First, order the points of the decorated dual lattice 
G in a manner consistent with the nested dissec- 
tion. The elements of the first dual separator A 
are at the end of this ordering. 

2. Using this ordering, set the values of the Kasteleyn 
matrix K , which is stored as a sparse matrix. 
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3. Factorize K using Pfaffian elimination. We use a 
pre-computed list of elementary operations to carry 
out the cross operations for all elements that are po- 
tentially non-zero. (Stop here if only the partition 
function for U is needed; the partition function Z 
is just the product of alternate superdiagonal ele- 
ments in Y, i.e., Z = Hfllo^ U^-) 

4. Using this factorization, compute the elements of 

that are indexed by elements of A; this is 
[K~^A, the lower-right submatrix of K^^ with in- 
dices contained in A. (For some speed up, as sug- 
gested in Ref. [13, we only compute the elements of 
[A'^i]^ that are needed in the following steps, at 
the time those elements are required.) 

5. Assign dimers (gi, ri), (q2, ^'2), • ■ • along the separa- 
tor A: 

(a) Choose a node gi e A such all edges 
that are incident upon qi arc fully con- 
tained in A. Choose among the po- 
tential edges (qi,ri) = e with the 
probabilities if(gi, ri) Pf (iC,, Pf(is:) = 

qi,riJ- 

(b) Repeat this last substep, [5al proceeding along 
the dimer and updating [K~^]q-^^ri,...,qk,rk ^nd 
its factorization, until no more matchings can 
be added wholly within the set A. 

6. Use the dimer matching for A to assign spin values 
in the spin separator D, which surroimd the dimer 
separator A. 

7. Recursively repeat items 1-7 for the subproblems 
U' and U". 

Note that, in some cases, an alteration of this proce- 
dure can be used to speed up this method. It might be 
that faster results can be obtained using simple float- 
ing point numbers (double precision), rather than multi- 
precision numbers, though they may not provide numer- 
ical accuracy to carry out all of the calculation. A com- 
promise would be to carry out the computation for only 
part of the separator at a time, making the computation 
more stable. The whole matrix K with the remaining un- 
chosen nodes is recomputed and the process is repeated. 
This method is asymptotically slower, but practical for 
systems of intermediate size at intermediate temperature. 

A. Entries of K: nested dissection and storage 

The Kasteleyn matrix K, as defined in Section III Al 
is indexed by the nodes of the decorated dual graph G. 
As the entries K{q, r) = ±a;(g, r) of K are non-zero only 
for entries indexed by neighboring points q and r on the 
decorated dual lattice, this 0{n) x 0{n) matrix has only 
0{n) non-zero entries. If the nodes are indexed in a nat- 
ural, geometric, lattice order, the Kasteleyn matrix K is 



simple, as shown in Fig. lUJb). However, matrix manipu- 
lations, such as Pfaffian elimination, for general matrices 
might lead to the computation of O(n^) non-zero entries. 

To compute the correlations between spins on the sep- 
arator, the nodes are reordered, though kept together in 
city groups. In this reordering, the nodes are each as- 
signed a new index. This reordering satisfies the nested 
dissection property that, at each level, the separator 
nodes in A, which give the spin sub-sample U, have the 
highest index. This implies that the non-zero values de- 
fined by the weights contained within the separator A are 
at the lower and rightmost parts of i^, at each level, while 
the non-zero values for nodes belonging to U' and U" [see 
Eq. ([7])] arc confined to square blocks about the diagonal. 
An example of the distribution of matrix entries, given 
this ordering of the nodes V of G, is shown in Fig. IDJc). 
This organization confines all matrix manipulation to a 
portion of the shaded regions of the matrix and to a nar- 
row band around the diagonal, as unshaded entries away 
from the diagonal always have value zero. The shaded 
regions make up 0{N'^/'^) entries, though only a subset 
of even those entries, growing with N approximately as 
^ iV, possibly with a logarithmic correction, arc used in 
the Pfaffian elimination. 

Given our specific choice of separator, the nodes of G 
corresponding to the Kasteleyn cities always form sub- 
sequences in the ordering of the nodes. That is, they 
remain grouped together. Note that the submatrices for 
each city are uniform in structure. This choice of sep- 
arator A (as all of the dual nodes between two rows or 
columns of spins) is not the most efficient, as slightly 
smaller separators A C G can be chosen, but it is a very 
convenient choice that maintains a uniform structure. 

We use this ordering to construct K as a sparse ma- 
trix, using 0{N) operations and time. The sparse ma- 
trix storage scheme is relatively direct (see, e.g. [i^ for a 
discussion on sparse matrix algorithms and storage tech- 
niques). We have the advantage here that, for the Pfaf- 
fian elimination, both the locations of the needed ele- 
ments and the list of operations using these elements can 
be pre-computed and stored on disk. This allows us to 
place the elements of the matrix K in a linear array with 
0{N) elements, with the elements ordered by the step 
at which they are first needed in the Pfaffian elimina- 
tion. This precomputation is independent of both the 
data type that we use and the bond strengths for the 
spin lattice. 



B. Pfaffian factorization 

Pfaffian elimination and the concomitant factorization 
of K proceed by the elimination of elements by cross op- 
erations. There are two types of cross operations that 
are carried out. The first type of operation eliminates 
all but the first of the non-zero entries in an even row. 
This is done for an even row q by using (see Eq. [5]) 
a{q,r) = -K{q,r)/K{q,q + 1) for all r > q + 2. The 
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Figure 5: Example of cross operations used for PfafHan elim- 
ination, (i) A skew-symmetric matrix K. (ii) The result of a 
cross operation K — > L{a,i, j)KL^ {a, q,r) of the first type, 
with q — 0, r = 2, a = —a, applied to K. This is found 
by adding a — —K{q,r)/K{q,q + 1) times column q + 1 to 
column r and then a times row q + 1 to row r, to eliminate 
the element at location {q,r). (iii) The result of the next 
cross operation, with g = 0, r = 3, and a = —b. (iv) The 
result of two subsequent operations of the second type, where 
a{q,r) — —K{q,r)/K{q,q— 1), for q = 1, r = 2 and q = 1, 
r = 3. These latter types of operation are not needed to com- 
pute Pf(-ft'), but are needed for finding [-K^'^Ja- The Pfaffian 
of K is the product of the superdiagonal elements in even 
rows: here, Pi{K) = (1)(1 -ac+b). 



second type eliminates all entries in odd rows. This is 
done for odd i using a{q, r) = —K{q, r)/K{q — 1, q) with 
r > q. Examples of operations of each type are traced 
out in Fig. [S] 

We note that in carrying out Pfaffian eUmination, a 
potential danger would be that one of the even-row su- 
perdiagonal elements, K{q,q + 1) with q even, is zero. 
In this situation, it would be necessary to do a pivot- 
ing operation, which would destroy the nested dissection. 
However, given that the Kasteleyn cities remain grouped 
together, the sequential pairing of nodes (0, 1), (2, 3), . . . 
is always a matching. Hence the Pfaffian of any upper 
left portion of the Kasteleyn matrix, as we have arranged 
it, is non-zero, as the Pfaffian counts matchings (in a 
weighted fashion) , and there is always a matching for the 
upper left portion of the matrix of unit weight. This 
implies that all superdiagonal elements in the even rows 
must be non-zero. This provides a "built-in" version of 
the permutation of nodes to accommodate a matching 
that is given in Wilson's paper [l7j . In the periodic case 
(Sec. IIVP , for certain boundary weight choices at = 
(T = oo), when the bond strengths have uniform magni- 
tude, there can be "accidental" cancellations which will 
cause this procedure to fail, as the signed weight of a sub- 
matching can be exactly zero, even though the Pfaffian 



is non-zero. In this case, permutation of the remaining 
elements of the matrix (i.e., "pivoting") would be needed 
to remove a zero from the superdiagonal and obtain the 
correct factorization. 

The factorization found by Pfafhan elimination. Eq. 
(fT^ . then allows for the easy computation of the parti- 
tion function for the given sample, at the temperature 
used to set the elements of if desired. The Pfaffian 
of the original Kasteleyn matrix is simply the product of 
the even superdiagonal elements of 



7V-1 



Vi{K) ^Wvi 



(19) 



1=0 



Note that this is the procedure, computation of the Pfaf- 
fian of K using nested dissection, was used by Galluc- 
cio et al. (2^ 1 to compute the partition hmction. In 
that work, to compute the partition function at a given 
temperature, the arithmetic is carried out modulo prime 
integers, for a selection of prime integers. The partition 
function at that temperature is then reconstructed by ap- 
plication of the Chinese remainder theorem. The whole 
partition function as a function of (3 can be found by 
polynomial interpolation in exp(— /?). This full calcula- 
tion works only if the couplings Jij are restricted to small 
integer values, typically Jy = ±1. 



C. Sampling: inductively factorizing K ^ 

At this point, though one has the partition function 
(from the even superdiagonal elements of Y), sampling 
spin configurations requires a bit more work. The sam- 
pling can be carried out by using only the lower right 
hand corner [^"^Ja oi . This part of the matrix en- 
codes all the correlations between the spins in Z?, on the 
separator of the sample, via the correlations of dimer cov- 
erings of A. These correlations are used to make dimer 
(and then spin) assignments along the geometric separa- 
tor. The description in this subsection is based upon Wil- 
son's description and notation p7| . only with a change 
in the factorization method (Pfaffian vs. Gaussian). 

To assign a dimer covering inside the separator A, the 
algorithm proceeds through each of the edges in G that 
are wholly contained within the node set A and computes 
the probability that that edge is covered by a dimer, con- 
ditioned on earlier assignments of dimcrs in the separa- 
tor. The algorithm proceeds inductively by calculating 
the probabilities for placing the (/c-l- l)'st dimer using the 
results of the calculations for the previous k edges in A, 
P = {{qi,ri),...,{qk,rk)}. 

The inductive computation of the probabilities are 
based on Eq. (jl7p . which in turn requires the compu- 
pf (i^ 

tation of the ratio ^pjpf^- This ratio is found from 

the change in the Pfaffian of [K~^]q^ ri,...,qk,rk that re- 
sults from the augmentation of [/'iT"^] by two rows and 
columns, those with indices qu+i and rk+i in . To 
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calculate this change, the algorithm maintains a factor- 
ization of = [/C~^]p which is tentatively updated to 
test the addition of an edge. This factorization allows for 
the ratios of PfafFians to be quickly computed. The ma- 
trix [-/^^^Ja is first found by computing a subset of the 
rows and columns of K^^ using Eq. (fT4|l and the Pfaffian 
factorization of K, Eq. (fT^ . 

To select matched edges within A, one considers in 
turn nodes q €z A such that all neighbors r of q are also 
in A and selects one of these neighbors with the correct 
probability. When considering matches for such a node 
Qk+i, assume that one has already selected k dimers in 
A, as part of a sampling inside A, and that one knows 
the matrices Mk and Vk in the factorization 



MkAuMl = Vk, 



(20) 



where all matrices in this equation are of dimension 
2k X 2fc, Mk is lower triangular, and Vk has the same 
super-diagonal structure as Y . For a given trial edge 
(qfc+i,rfe+i), we can tentatively extend the matrices Mk 
and Vk to Mk+i and Vfc+i, with 



'J-k+l 



and 



14+1 



Mk 
mfc+i / 



Vk 

Vk+1 



where -I- 1 is a 2 x 2 antisymmetric matrix. 



Vk 



Zk+i 
^Zk+i 



(21) 



(22) 



(23) 



and TUk+i is a 2 X 2fc matrix. To compute these trial 
solutions Mk+i and V^+i, one first tentatively updates 



Ak+i = 



Ak 



(24) 



using the rows and columns indexed by qk+i and rk+i 
from [if~"'^]A to fill in Ak+i and reading off ak+i and 
hk+i ■ Direct matrix multiplication and requiring Eq. (|20p 
for Ak+i then give 



TTik+l 

and that 



-ak+iA], 



-ak+iM^V^H'Ik 



Zk+l 



= bk 



k+l 



ak+iMlV-HlkO^ 



k+l 



(25) 



(26) 



As Pf(Afc) — Hi^i fc^ii the factor Zk+i is the ratio 
'Pi{Ak+i) / ^i{Ak) of the Pfaffians that is needed to apply 
Eq. ([T7)) . Hence, this update in the factorization allows 
us to find the probability Xq^^^^rk+iZk+i{qk+i,rk+i) of 
selecting the specific edge {qk+i,rk+i) to augment the 
matching. Once we have chosen a match for qk+i, we 
then update Ak to Ak+i from , Vk using z^+i, and 



k=0 



k=l 



k=2 



k=3 



k=3 



k=4 



Figure 6: [color online] An example of the dimer assignment 
procedure for a separator A that is three cities wide. Initially, 
no edges are matched (left part of = 0). The first choice, 
A; = 0, is between the two edges inside A that are incident 
upon the far left node. In the example shown, the lower bond 
(connecting node to node 2; see Fig. [Sjl is chosen, as shown 
on the left of the fc = 1 section of the figure. At this stage, 
one has computed matrices A\, Mi, and Vi. The comparison 
of the next two possible matchings, shown on the right part 
of the k = 1 subfigure, compares the inclusion of the (g, r) = 
(3, 5) and (3, 4) edges. In some cases, as in the first k = 3 
panel, a choice is forced and Ak, Mk, and Vk need not be 
updated. The = 4 choices are forced, as an even number 
of domain walls must cross the separator, so that an even 
number of the top nodes and an even number of the bottom 
nodes are unmatched. 



Mk using Eq. (^5]) . This process is repeated until a max- 
imal (though usually not complete) matching within A 
is obtained. With our choice of Fisher cities, there are 
only two candidates Vk+i for each qk+i when using fixed 
boundary conditions; for periodic boundary conditions 
(Sec. IIV[) . matching the initial node qi = requires the 



comparison of three choices. Note that not all the z need 
be computed as the total probability sums to unity; when 
considering two choices, considerable time is saved by 
computing the probability of only one of the choices. An 
example of dimer assignment is depicted in Fig. [6l 

The results derived by Wilson for the bounds on the 
number of steps using Gaussian elimination carry over 
directly to the approach using Pfaffian elimination. The 
maximal size of the separator is of order 0(V) = 0{n^^^). 
There are 0(n'^/^) operations in the dimer assignment for 
the largest separator: matching a single dimer requires 
at most 0{n) steps, due to the multiplication of matrices 
of size 2 X 0{n^^^) by matrices of size 0{n^/^) x 0{n^/^), 
and there are 0(n^/^) matchings in each separator. Cal- 
culating K^^ also requires 0{n^^^) steps. As the smaller 
separators decrease in size geometrically, as the sample 
is subdivided, the number of operations for each of the 
smaller separators decreases geometrically, and the sum 
of steps over all levels of the nested dissection gives a 
total of 0{n^^^) arithmetic steps to generate a random 
assignment. The running time then is a product of the 
time per operation, which depends on the needed pre- 
cision, and this number of steps. As discussed in more 
detail in Sec. HITEI the running time grows roughly lin- 
early with the precision: the necessary precision grows 
only slowly with n, but proportionally to f3. 

Once all nodes in the separator A have dimers asso- 
ciated with them, the broken bonds along the strip D 
of the Ising system are found from the locations of the 
dimers between these cities and the neighboring ones. 
We can then directly assign the spins along the strip. An 
example of such a spin assignment is displayed in Fig. [T] 



D. Verification 
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(a) (b) 




Figure 7: [color online] Application of the result of the exam- 
ple dimer assignment from Fig. [6] to the spin assignment, 
(a) The initial spin configuration, with fixed spins on the 
boundary. The portion of G used to compute K, [K~^]a 
and the dimer assignment is indicated in gray. The middle 
row of Fisher cities composes A, the dimer separator, (b) 
The sample dimer assignment (partial matching) for A from 
Fig. |6] superimposed on G. (c) Extra choices in the matching 
are forced by the matching internal to A. These additional 
dimers cut across the bonds separating spins in D, the two 
spin rows parallel to A. (d) In the last step, the modifiable 
spins are updated. The update is based upon the portion of 
domain walls forced by the partial matching in (c). Moving 
from left to right, for example, from the two fixed spins on 
the middle of the left side, a spin is reversed if an odd number 
of dimers extending from A are crossed. 



The structure of the calculation is rather complex, so 
we verified our implementation of the algorithm in several 
ways. We checked exact partition function calculations 
for pure systems against the results of our computation. 
Exact enumeration for pure and disordered samples in 
systems up to n = 5^ was used to predict sampling prob- 
abihties: we then used our code to generate over 10^ 
samples and compared the sampled probability distribu- 
tion with the exact calculations. These were in statistical 
agreement. Each author of this paper developed a code 
independently: these were compared on the same Gaus- 
sian spin glass samples of size 33^ and found to generate 
the same distribution for configurations, at low tempera- 
tures, also consistent with the Boltzmann distribution for 
total energy. At low temperatures, the sampled configu- 
rations approached those of the ground state configura- 
tions (which were predicted using an independent ground 
state code based on combinatorial optimization methods 



E. Data types and timing 

Our code is constructed so that the data type of ma- 
trix elements can be any field (double precision numbers, 
multiple-precision numbers, or exact rationals, for exam- 
ple). This allows us to check the effects of the choice 
of numerical type on the accuracy, stability, and run- 
ning time of the sampling algorithm. For higher pre- 
cision variables, we use the GMP library 14611 for exact 
rational arithmetic and cither the MPFR |43| extension 
to GMP or the GMP library itself for multiple-precision 
floating point arithmetic. We find that the latter two 
floating point types give comparable performance and ac- 
curacy. Using exact rationals allows for mathematically 
exact sampling, but results in a temperature-dependent 
slowdown by a factor of 10 or 100 over the range of tem- 
peratures, T = 0.1 to T = 1, we used while comparing 
rationals with floating point calculations. 
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Figure 8: [color online] The plot shows the sample averages 
of the number of bits of precision Baiig required to obtain the 
correct sampling, for system sizes n = = 6^ through 126^ 
with fixed boundaries, as a function of inverse temperature 
/3, for Gaussian disorder. The lines indicate linear fits of the 
form Ba_vg = c(5. The number of bits needed for periodic 
boundary conditions (not shown) are very close to these same 
lines, at each system size. To find an accurate result with 
high confidence, one can use twice the average needed value; 
this was sufficient for all samples (> 10^) that we examined. 



The edges u are chosen by comparing the probability 
P{p,u \p) with a random number chosen in the interval 
[0, 1). The sequence of random numbers and computed 
probabilities determines the spin configuration selected. 
We determine the needed precision for a given sample 
and temperature by demanding that the result of a spe- 
cific assignment be independent of the precision, for a 
given sequence of random numbers. Note that using this 
precision does not give the exact values of the probabil- 
ities at each stage of the computation, but the sampling 
does not change at increased precision. If a number in 
the sequence happens to be extremely close to the com- 
puted probability, higher precision arithmetic could be 
required. 

Results of our tests for needed precisions arc summa- 
rized in Fig. [51 where we plot the number of bits needed, 
determined by bisection in the number of bits, averaged 
over random number sequences and disorder. We find 
that the distribution of the required number of bits is 
not very broad, regardless of temperature and disorder 
realization J'. Less than 10""* of the attempts require 
more than double the average precision to find the cor- 
rect sampling. Hence fixing the precision at two to three 
times the average value will almost guarantee an exact 
sampling. 

For high temperatures (of order T = 1), low precisions 
(i.e. fixed double precision variables) are sufficient for the 
system sizes we study (see Fig. For lower tempera- 
tures, higher precisions are needed. The needed precision 
is well fit by a linear growth in /3, for f3 > 0.5. This is 
consistent with the expectation that, as the weights vary 
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Figure 9: [color online] Run time, measured in seconds, to 
generate a single configuration, as a function of system size L, 
using a 2.4 GHz Intel Core 2 Duo processor (MacBook Pro). 
Double precision (64 bit floating point) data is indicated with 
triangles, while multi-precision results for B = 512, 2048, and 
8192 bits are indicated by squares, diamonds, and circles, re- 
spectively. Samples are generated for L < 128 with fixed 
boundaries (closed symbols) and for L < 64 with periodic 
boundaries (filled symbols). The sample-to-sample fluctua- 
tion of disorder realization is less than 0.1% of the run time, 
so error bars are not shown. The solid line indicates the form 
of the expected dependence of run time on system size for a 
given flxed precision, that is, ~ L'^ . 



as exp(— /3J), the number of bits needed to describe the 
weights grows linearly with /?, for fixed typical values of 
J. The number of needed bits grows only slowly with L. 
This is consistent with the structure of the sampling and 
Pfaffian computation, which are hierarchical in structure, 
so that the accumulated error grows only slowly with L. 

For systems up to size 64^, 600 bits of precision are 
sufficient for temperatures T > 0.1. For larger systems 
and lower temperatures, more bits are needed. For ex- 
ample, we use 2048 bits to reliably sample configurations 
at (3^25 and L = 128. 

We collected timing data for the performance of our 
algorithm as a function of system size and temperature. 
These data are summarized in Fig. [9l We find that sam- 
pling with periodic boundary conditions (Sec. IIVP takes 
approximately 5.5-6.5 times longer than sampling with 
fixed boundary conditions. The needed precision and 
running times for ±J disorder arc very close to those 
shown in Figs. [9] and [H For 64 < B < 512, the run time 
to sample a configuration varies only slowly with B, ap- 
proximately by a factor of 1.5 over this range. For higher 
precision, the running time grows somewhat faster than 
linearly with B, and hence somewhat faster than linearly 
with f3. 



IV. PERIODIC BOUNDARY CONDITIONS 

Fixed boundary conditions arc appropriate for patch- 
work dynamics, but, for other simulations, other bound- 
ary conditions, may be useful. One simple way to imple- 
ment open boundary conditions is to set to zero all the 
Jij connecting interior to boundary spins. For cylindri- 
cal samples with open boundaries, we use a "separator" 
which does not actually separate the graph, but one that 
slices the sample perpendicular to the circumfcrcncc of 
the cylinder, resulting in a simple planar graph with fixed 
boundaries. Toroidal graphs require a more complicated 
sampling scheme, as they are not planar. In general, for 
a graph of genus g, the partition function of the dimer 
problem may be calculated exactly by summing 4^ Pfaf- 
fians [48l |. The reasoning behind this summation can be 
adapted to sampling for periodic spin lattices. 



A. Partition function on the periodic lattice 

The Kasteleyn matrix approach for computing Z can 
be extended to handle the periodic case, by adding con- 
nections between cities that complete the periodic bound- 
aries, converting the planar square sample to a toroidal 
one, but the direct correspondence between dimer con- 
figurations and spin configurations is affected. On the 
torus, topologically non-trivial domain walls must always 
come in pairs, or the spin configuration can not be con- 
sistently defined. But the matching problem allows for 
odd numbers of loops to wrap around the torus on cither 
axis. For T = ground states, one can decide to ignore 
this fact and allow variable boundary conditions, which 
allow for an odd number of domain walls relative to other 
boundary conditions. Choosing the boundary condition 
and spin configuration that jointly minimize Ti. gives the 
extended ground state construction [l^. At finite tem- 
perature with fixed boundary conditions, however, we 
need to arrange for the cancellation of dimer configura- 
tions which would imply an odd number of domain walls 
that wrap around either axis. 

This cancellation is achieved by summing over four 
Pfaffians, in a fashion similar to that developed for the 
primal lattice [l^, though the details differ for the dual 
lattice. The four Pfaffians correspond to four possible 
choices of sign for the elements of K{q,r) that complete 
the periodic connections. That is, the values of K{q, r) 
for edges that connect the last column to the first col- 
umn (that wrap around in the x direction) arc uniformly 
set to one of two choices, ± exp[— /3w(q, r)], and the val- 
ues for the edges that connect the last row to the first 
row (that wrap around in the y direction) are also uni- 
formly set, independent of the choice for the cc-wrapping 
bonds, again to ± exp[— /3w((7, r)]. This gives four matri- 
ces, K'^'^ , K K'^ , and K . The dimer configu- 
rations that arc summed up in the Pfaffians enter with 
different relative signs, depending on how many times 
the matchings wrap around each axis, as the parity of the 
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Table L A table of the signs for different combinations of 
spanning loop parities in the dimer model for each of the four 
Pfaffians for the torus. The set of loops found from a 

dimer configuration can have a total wrapping number that 
is odd (o) or even (e) number along either the horizontal or 
vertical directions. This gives four possible classes of dimer 
configurations (e,e), (e,o), (o,e) and (o,o). For the dual map- 
ping used here, the physical spin configurations for the Ising 
model are restricted to those with an even number of domain 
walls wrapping in both directions, i.e., the (e,e) class. The 
four classes of dimer configurations are summed in each Pfaf- 
fian of the four Kasteleyn matrices, , with a sign that de- 
pends on the class and the matrix. These four matrices assign 
difi'erent signs to the weights of the dual edges that connect 
the boundaries together, with a -I- or — sign for each of the 
two types of boundary connections, i.e., horizontal or vertical. 
Applying this table, we get the partition function for the valid 
dimer configurations by the sum Z = [Pf(A++)-hPf(A+-)-H 
Pf(A"+) -I- Pf(A )]/2, which counts only the (e,e) class 
of dimer configurations. This sum differs from the more com- 
monly studied case, the dimer model using cities on the primal 
lattice, where all classes of matchings are valid configurations 
and Z = [- Pf(A++) + Pf(A+-) + Pf (A-+) + Pf(A— )]/2 
gives the sum over (e,e), (o,e), (e,o), and (e,e). 

windings affects the sign of the dimer configurations when 
the negative sign is chosen for the periodically-connecting 
edges. The effects of these signs are tabulated and ex- 
plained in Table [J The sum of the Pfaffian of these four 
matrices then gives twice the partition function, as those 
dimer configurations with an even number of wrapping 
loops enter four times and those with an odd number, 
in either direction, are cancelled out, and there is a two- 
to-one mapping of spin configurations to domain walls in 
the periodic case (due to global spin flip symmetry). 



B. Matching probabilities for the torus 

There are several simple possible choices for a nested 
dissection for toroidal samples of dimension L x L. The 
number of cities is the same as the number of variable 
spins, i.e., L^. We chose to use a horizontal strip of length 
L in the first row of cities, which fixes the spins in the first 
two rows, followed by a vertical strip in of length L — 1 
in the first column, which fixes the first two columns of 
spins, followed by a sampling the remaining (i— 1) x (L — 

1) cities, i.e., a sampling of the remaining (i — 2) x (i — 

2) spins using the already determined spins in the first 
two columns and rows as fixed spin boundary conditions. 
The first two "separators" don't divide the sample into 
separate pieces, but instead provide for the cutting of 
loops that wind around the torus, in two stages. 
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For the first sampling, the periodic horizontal row, one 
has to sample using four Kasteleyn matrices in parallel. 
For the second sampling, on a cylindrical geometry, one 
needs to find probabilities by summing over two Kaste- 
leyn matrices and ^ in order to eliminate domain 
walls that wrap around the cylinder an odd number of 
times. We can consider both cases as specific examples 
of a general problem: sampling using multiple Kasteleyn 
matrices simultaneously. 

For this general case, consider a partition function Z 
that is found by summing the PfafBan over matrices A"", 
with weights qa (e.g., a = ±± and qa ^ \ for toroidal 
boundary conditions). The partition function is then 

a 

The computation of probability of selection is more com- 
plicated than for the case of a single K. For each AT", 
we consider the inverse indexed by elements of the sepa- 
rator A, [(A'")~^]a, and inductively factorize [{K°')~-^]p 
for our current choice of sampled edges p = {ei, . . . , e^-}. 
The conditional probability of choosing edge Cfc+i, sim- 
plifying the notation by writing u for Ck+i and using 
2;"(e) to denote for edge e = (gfe, r^), is then given by 

Pip,u\p) = ^ (28) 

Zjp 

^ E„gaPf(A7,jnegp..^"(e) 

E„9.Pf(i^p")neep^"(e) 

^ Eg la miK'^rXu) Pf(A-) Ueep,u ^"(e) 

E„ P/([(K")-i]p) Pf(i^") Deep ^He) 
^ EagcPf(^")nesp,.^°(eK(e) 
E„9„Pf(i^")ne6p^"(e)2:"(e) 

^ EaC«(p)^"(^)a:°(u) 
EaCa(p) 

where 

Up)=mKnnJ-^l^^. (29) 

This extra weighting quantity, Ca{p), is not needed for 
planar samples, due to cancellations, but is required here 
to allow for the different p-dependent weightings resulting 
from the distinct boimdary conditions. It incorporates 
the weight of the whole A"" matrix, the modification of 
those weights by the factors of z"(e) resulting from the 
choice of edges in p, and the sign of the weights (the mag- 
nitudes are identical in each a for a given choice of p and 
hence cancel out). This weighting factor can be updated 
at each stage k along with the set of , , and A'^ for 
each a. In the case of the periodic lattice, these four sets 
of matrices are updated and used to compute the values 
of z"(e) to find the conditional probabilities. 




Figure 10: Relative domain walls found in an individual 2D 
Ising spin glass sample with 64^ spins, periodic boundary con- 
ditions, and unit variance Gaussian disorder, for temperature 
T = 0.16. The bond satisfaction probabilities were estimated 
by averaging over 660 samples. As in Fig. [TJ the lines indicate 
the probability of relative domain walls between two config- 
urations: the darkest lines indicate where the bond dual to 
that domain wall has a nearly 50% chance of opposite or equal 
relative orientations; where there is no line or a light line sep- 
arating two spins, the two spins have a very high probability 
of a single relative orientation, either aligned or opposite. 

C. Sampling spins 

The dimer assignments are carried out on G for the 
periodic case using Eq. ((28|) . To finally carry out the 
sampling on the torus, one first arbitrarily sets the value 
of an initial spin, the spin at the upper left corner, i.e., 
at location (0, 0). The spin at the left side of the second 
row, at location (1,0), is fixed by the first element of the 
matching for the first separator. This is the exceptional 
case for this lattice where one has three choices for the 
matching edge on G ((0,1), (0,2), and (0,6L-2)). After 
this choice has been made, the rest of the spins in the first 
two rows are then assigned as in the fixed boundary case. 
An example of the relative domain wall density for a 64^ 
periodic sample is displayed in Fig. 1101 This plot shows 
the variance /iy (1 — /iy ) in the bond satisfaction, where 
fiij is the probability of a given bond being satisfied, i.e., 
siSjJij ^ 0. 

D. Running time 

We find that the number of bits required for the pe- 
riodic case increases only by a small amount, about 1%, 
over the planar case for samples of the same size. Car- 
rying out the initial Pfaffian elimination for single a for 
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the entire sample is slower than for the planar case, as 
there are about four times as many operations, but this 
computation requires only a small fraction of the time 
in any case. However, as the periodic case requires the 
maintenance of four Vk, Ak, and Mk matrices, sampling 
in the periodic case is slower than for the fixed bound- 
ary case. We find that sample generation is about 5.5 
times slower for periodic samples, compared with planar 
samples, for L = 16 through L = 64. 

V. CONCLUDING COMMENTS 

In this paper, we have described an algorithm that 
generates spin configurations for the 2D Ising spin glass, 
where the samples generated are directly selected ac- 
cording to the equilibrium probability distribution. This 
method follows from Wilson's dimer sampling algorithm, 
though we have modified the matrix algebra for speed 
and simplicity, and have adopted the dimer matching to 
the study of the Ising spin glass. We have also generalized 
the method to periodic samples. 

We note that as the inverse Kasteleyn matrix contains 
the dimer-dimer correlation functions along the separa- 
tor, one need not carry out all of the sampling steps to 
compute domain wall densities. One can directly ex- 
amine the inverse on the separator to find the domain 
wall densities on a single separator, by stopping at step 
|4] of the outline in Sec. IIIII The separator can then be 



changed to compute the bond satisfaction probabilities in 
each row of the sample. Sampling configurations provides 
more information, but if the bond satisfaction variance 
is all that is needed, this approach is more precise and is 
not unreasonably slow. 

This algorithm can also be used to directly and uni- 
formly sample ground states in the 2D ±J spin glass 
model. At low enough temperatures (on the order of 
T w 0.1), the ground states occur frequently, as can be 
confirmed by their energies being lowest or by compari- 
son with a ground state energy found by combinatorial 
optimization. The statistics of the ground state config- 
urations can therefore be directly sampled (by rejecting 
other states when they occur), exactly, using this algo- 
rithm. 

Our implementation of the sampling algorithm is effi- 
cient enough to allow for rapid enough sampling to study 
finite temperature patchwork dynamics out to patch sizes 
£ of at least £ = 32. Large numbers of samples can be 
comfortably generated for L = 64 and T < 0.01 or for 
L = 128 and T = 0.04. This should allow for more 
conclusive studies on the Gaussian and ±J spin glass 
problems in two dimensions. 

We thank Jan Vondrak for sharing his code for com- 
puting the partition function of the ± J spin glass model 
and Simon Catterall for sharing his code for computing 
Pfaffians. This work was supported in part by the NSF 
under grant No. 0606424. 
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